CAM Implementation of Mesoscale eddy diffusivity based on mixing theory
Here are my notes on the development of a parameterization for eddy diffusivity. The scheme is based on the paper by Michael Bates, Ross Tulloch, John Marshall, AND Raffaele Ferrari. I'm working with Gokhan Danabasoglu and Matt Long here at NCAR and John Marshall at MIT.
Gokhan's Implementation Notes on Diffusivity with Steering Level Suppression
12 November 2014
Starting with $$K=u_{rms}∗L_{mix}$$
the general form for diffusivity, $K$, is given by equation (6) of [Bates et al. (2014)] http://journals.ametsoc.org/doi/abs/10.1175/JPO-D-13-0130.1. Namely,
$$L_{mix} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $$
Here,
$u_{rms}$ is the root-mean-square (rms) eddy velocity;
$L_{mix}$ is the mixing length;
$L_{eddy}$ is the eddy diameter (depth independent);
$u_{mean}$ is the mean zonal velocity (resolved);
$c$ is the zonal eddy phase speed (depth independent);
$\Gamma = 0.35$;
$b1 \sim 4$.
Eddy Length Scales:
Rossby deformation radius = $L_r = {c_r \over |f|}$
Equatorial Rossby deformation radius $= L_{req} = {\sqrt{c_r \over 2\beta}}$
Rhines scale $= L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma_{vi} \over \beta}$
$c_r$ is the first baroclinic wave speed computed following equation (2.2) of Chelton et al. (1998) with $m=1$;
$f$ is the Coriolis parameter;
$\beta$ is the latitudinal variation of the Coriolis parameter; and
$\sigma_{vi}$ is the Eady growth rate given by
$\sigma_{vi} = {f \over \sqrt{R_i}}$
with $R_i$ the vertically integrated (over the 100 – 2000 m depth range) Richardson number.
So, any one these length scales could be used as an eddy length scale. An alternative is
$L_{eddy} = min (L_r, L_{req}, L_{Rh}).$
Eddy Velocity:
$u_{rms} = alpha*\sigma_{vi}*L_r$
where $\sigma_{vi}$ is the Eady growth rate based on local Richardson number and $\alpha$ is a scaling constant.
Zonal phase speed:
$c = - \beta * L_r^2$
References:
Bates, M., R. Tulloch, J. Marshall, and R. Ferrari, 2014: Rationalizing the spatial distribution of mesoscale eddy diffusivity in terms of mixing length theory. J. Phys. Oceanogr., 44, 1523-1540, doi: 10.1175/JPO-D-13-0130.1.
Chelton, D. B., R. A. deSzoeke, M. G. Schlax, K. E. Naggar, and N. Siwertz, 1998: Geographical variability of the first baroclinic Rossby radius of deformation. J. Phys. Oceanogr., 28, 433-460.
Tullock, R., J. Marshall, and K. S. Smith, 2009: Interpretation of the propagation of surface altimetric observations in terms of planetary waves and geostropic turbulence. J. Geophys. Res., 114, C02005, doi: 10.1029/2008JC005055.
Questions:
In the 2D implementation, $u_{rms}$ and $u_{mean}$ specifications: upper-ocean vertically or integrated or at $z = 0$?
A: $U_{rms}$ is not depth dependent; both $u_{rms}$ and $u_{mean}$ are for surface only
In the 2D implementation, vertical profile will be specified by $N2(z)$?
A: Yes, $N^2(z) \over N_{ref}(z)$ to be more precise
Local $R_i$ use imbedded in sigma in $u_{rms}$ calculation?
A: No, $u_{rms}$ is depth independent
alpha = ?
A: trial and error
Cancellation of $f$’s in $u_{rms}$ calculation?
Zonal phase speed equation correct? Both $\beta$ and $L_r$ will be positive, producing $c < 0$ always. This appears to be in contrast with Tullock et al. (2009).
A: What is plotted in Bates is (U-c)
Parameterizing the Eddy Length Scale
$$K=u_{rms}∗{\Gamma * \color{red}{L_{eddy}} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $$
The three length scales under consideration:
- Rossby Deformation Radius $\quad\quad L_r = {c_r \over |f|}$
- Rossby Equatorial Radius $ \quad\quad= L_{req} = {\sqrt{c_r \over 2\beta}}$
- Rhine's Scale $ \quad\quad L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma_{vi} \over \beta}$
The Rossby Deformation Radius $L_r$ depends on the the Baroclinic wave speed $c_r$ and the Coriolis force $f$. The CESM value of the first Baroclinic wave speed is derived as per eq 2.2 in Chelton 1998.
First Baroclinic Wave Speed $c_r$ Chelton vs CESM
# <!-- collapse=True -->
%pylab inline
%run steering_setup ;
basedir='/scratch/jet'
###nc = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nday1.yr.0249.nc')
nc = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nday1.0249-01.nc')
omega=7.2921e-5
deg2rad=pi/180
f=2*omega*sin(lat*deg2rad)
# <!-- collapse=True -->
### get Variable from file c_rossby is in cm/s convert to m/s
c_rossby = nc.variables['C_ROSSBY'][0,:,:]
#Resample (aka re-project, re-grid) the NCEP data to target grid. First with nearest neighbour resampling...
c_rossby_nearest = pyresample.kd_tree.resample_nearest(orig_def, c_rossby, \
targ_def, radius_of_influence=500000, fill_value=None)
c_rossby_mps = c_rossby_nearest/100.
fig1=MapZoneContour(lon_targcyc,lat_targcyc,c_rossby_mps,figsize=(14,9),
title="CESM First Baroclinic Wave Speed (m/s)",
fmt='%.1f' )
c_r_chelton=Image(filename=basedir+'/steering/chelton_sfig1.jpg')
display(c_r_chelton)
First Baroclinic Rossby Radius Chelton vs CESM
CESM Rossby deformation radius = $L_{eddy}$
$L_r = {c_r \over |f|}$
$L_{req} = {\sqrt{c_r \over 2\beta}}$
$\require{cancel} \cancel {L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma \over \beta}}$
$L_{eddy} = min (L_r, L_{req}, \cancel{L_{Rh}}).$
# <!-- collapse=True -->
fcor = nc.variables['FCOR'][0,:,:]
fcor_nearest = pyresample.kd_tree.resample_nearest(orig_def, fcor, \
targ_def, radius_of_influence=500000, fill_value=None)
absf=np.abs(fcor_nearest)
l_rossby_calc=c_rossby_nearest/absf
btp = nc.variables['BTP'][0,:,:]
btp_nearest = pyresample.kd_tree.resample_nearest(orig_def, btp, \
targ_def, radius_of_influence=500000, fill_value=None)
l_rossbyeq_calc=np.sqrt(c_rossby_nearest/(2*btp_nearest))
l_eddy_calc=np.minimum(l_rossby_calc,l_rossbyeq_calc)
### convert regridded radius from cm to km
cm2km=1.e-5
l_eddy_calc_km=l_eddy_calc*cm2km
clevs=[10,20,30,40,50,60,80,100,150,230]
fig2=MapZoneContour(lon_targcyc,lat_targcyc,l_eddy_calc_km,figsize=(14,9),
levels=clevs,
title="CESM First Baroclinic Rossby Radius (km) $L_{eddy} \sim min(l_r,l_{req})$",
fmt='%.1i' )
l_r_chelton=Image(filename=basedir+'/steering/chelton_fig2.jpg')
display(l_r_chelton)
Parameterizing the Zonal Eddy Phase Speed $\color{red}{c}$
$K=u_{rms}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - \color{red}{c}|^2 /u_{rms}^2 (z=0)} $
$\cancel{c = - \beta * {L_r^2}}$ $L_r$ too high at equator
$c = - \beta * L_{eddy}^2$
# <!-- collapse=True -->
#beta
clevs=arange(1.e-14,25.e-14,1e-14)
fig,ax,cbar =MapContourg(lon_targcyc,lat_targcyc,btp_nearest,
levels1=clevs,
# norm=norm,
addzonal=True,
figsize=(14,5),
title="Beta (1/cm*s)")
l_eddy_calc_sq=l_eddy_calc*l_eddy_calc
clevs=np.logspace(1,15,40)
norm=LogNorm()
fig,ax,cbar =MapContourg(lon_targcyc,lat_targcyc,l_eddy_calc_sq,
addzonal=True,
levels1=clevs,
norm=norm,
figsize=(14,5),
title="$L_{eddy}^2 (cm^2)$ log scale")
#cbar.ax.get_xaxis().get_major_formatter().labelOnlyBase = True
#cbar.ax.get_yaxis().set_major_formatter(plt.LogFormatter(10, labelOnlyBase=True))
#ax.yaxis.set_major_locator(ticker.MultipleLocator)
#ax.yaxis.set_major_formatter(LogFormatter())
#cbar.formatter.set_scientific(False)
#cbar.formatter.set_powerlimits((0, 0))
#cbar.update_ticks()
cbar.set_ticks(np.logspace(1,15,15))
c=btp_nearest*l_eddy_calc*l_eddy_calc # units=1/(cm s) * cm^2 = cm/s
clevs=arange(-10,20,1)
fig,ax,cbar = MapContourg(lon_targcyc,lat_targcyc,c,
addzonal=True,
levels1=clevs,
# norm=norm,
setover='darkred',
figsize=(14,5),
title="eddy phase speed c = beta*leddy**2 (cm/s) linear scale(NOTE I'm plotting positive c here - changed to negative after this plot")
cbar.set_ticks(np.linspace(-10,20,7))
c=-1.*c
Hughes Phase Speed (cm/s) from Tulloch Marshall Smith '09
# <!-- collapse=True -->
hughes_c=Image(filename=basedir+'/steering/Hughes_phase_speed.png')
display(hughes_c)
Parameterizing the Zonal Velocity $\color{red}{u_{mean}}$
$K=u_{rms}∗{\Gamma * L_{eddy} \over (1 + b1 * |\color{red}{u_{mean}} - c|^2 /u_{rms}^2 (z=0)} $
# <!-- collapse=True -->
################# PLOT (u-c)^2 convert m^2/s^2 #######################
umean=nc.variables['U_MEAN'][0,:,:]
umean10_nearest = pyresample.kd_tree.resample_nearest(orig_def, umean,
targ_def, radius_of_influence=500000, fill_value=None)
umean10_mps=umean10_nearest*1.e-2
clevs=np.linspace(-60,60,40)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,umean10_nearest,
addzonal=True,
figsize=(14,5),
levels1=clevs,
# setover='darkred',
# setunder='darkblue',
title="U zonal velocity (cm/s)")
#fig.subplots_adjust(left=0.15, right=.9, bottom=0.1, top=0.9);
cbar.set_ticks(arange(-60,60,10))
pylab.xlim([-20.,20.])
# get 2D versions of the lat and lon variables add longitude start point here!
lon_ecco_orig, lat_ecco_orig = np.meshgrid(lon_ecco[:], lat_ecco[:])
orig_ecco_def = pyresample.geometry.GridDefinition(lons=lon_ecco_orig, lats=lat_ecco_orig)
ecco_u_nearest = pyresample.kd_tree.resample_nearest(orig_ecco_def, ecco_u, \
targ_def, radius_of_influence=500000, fill_value=None)
ecco_u_nearest_masked=np.ma.masked_where(ecco_u_nearest<-10, ecco_u_nearest, copy=True)
#print ecco_u_cyc_masked[:,0]
clevs=np.linspace(-60,60,40)
#ecco_u_nearest = pyresample.kd_tree.resample_nearest(orig_def, ecco_u, \
# targ_def, radius_of_influence=500000, fill_value=None)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,ecco_u_nearest_masked*100,
addzonal=True,
figsize=(14,5),
levels1=clevs,
# extend='both',
title="ECCO U zonal velocity (cm/s)")
cbar.set_ticks(arange(-60,60,10))
pylab.xlim([-20.,20.])
# <!-- collapse=True -->
clevs=-np.logspace(3,-3,20)
norm=matplotlib.colors.SymLogNorm(linthresh=1e-4, vmin=-1000, vmax=-.001)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,c,
addzonal=True,
levels1=clevs,
norm=norm,
figsize=(14,5),
title="eddy phase speed c = beta*leddy**2 (cm/s) log scale")
# Customize y tick lables
yticks=[-1000,-500,-200,-100,-50,-20,-10,-5,-2,-1,-.500,-.200,-.1,-.05,-.02,-.01,-.005,-.002,-.001]
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
#ax.ticklabel_format(axis='y', style='sci', scilimits=(-2,2))
#majorFormatter = ticker.FormatStrFormatter('%.3f')
#majorLocator = ticker.LogLocator(base=10.0, subs=[1.0], numdecs=6, numticks=15)
#ticker.MultipleLocator(20)
#cbar.ax.yaxis.set_major_locator(majorLocator)
#cbar.ax.yaxis.set_major_formatter(majorFormatter)
#ax.yaxis.get_major_formatter().set_powerlimits((0, 1))
pylab.xlim([-200.,0.])
uminusc=umean10_nearest-c #units cm/s
clevs=np.linspace(0,25,25)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc,
addzonal=True,
# norm=norm,
levels1=clevs,
figsize=(14,5),
title="U-c (cm/s) (max = %i/min = %i) (c calculated using Leddy)"%(uminusc.max(),uminusc.min()))
yticks=linspace(0,25,6)
zonalticks=np.linspace(0,20,3)
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
pylab.xlim([0,25])
pylab.xticks(zonalticks)
# <!-- collapse=True -->
bates_uminusc=Image(filename=basedir+'/steering/bates2013_uminusc.png')
display(bates_uminusc)
# <!-- collapse=True -->
uminusc_mps_sq=uminusc*uminusc*1.e-4 #units cm/s*cm/s*m/100cm*m/100cm
clevs=np.logspace(-4,1,40)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc_mps_sq,
addzonal=True,
levels1=clevs,
figsize=(14,5),
norm=norm,
title=" ${(U-c)}^2 {(m^2/s^2)}$ max = %g/min = %g (log scale) "%(uminusc_mps_sq.max(),uminusc_mps_sq.min()))
#customize colorbar ticks
yticks=[1e-4,3e-4,1e-3,3e-3,1e-2,3e-2,1e-1,3e-1,1,3,10]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06])
clevs=np.linspace(0,.06,25)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc_mps_sq,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title=" ${(U-c)}^2 {(m^2/s^2)}$ max = %g/min = %g "%(uminusc_mps_sq.max(),uminusc_mps_sq.min()))
#customize colorbar ticks
yticks=[0,.02,.04,.06]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06])
# <!-- collapse=True -->
bates_uminuscsq=Image(filename=basedir+'/steering/Bates_u_minus_c_sq.jpg')
display(bates_uminuscsq)
Parameterizing the root-mean-square eddy velocity $\color{red}{u_{rms}}$ and the Eady Growth Rate $\color{red}{\sigma_{vi}}$
$K=\color{red}{u_{rms}}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $
$u_{rms} = alpha*\color{red}{\sigma_{vi}}*L_{eddy}$
The original derivation of $\sigma_{vi}$:
$\sigma_{vi} = {f \over \sqrt{R_i}}$
with $R_i$ the vertically integrated (over the 100 – 2000 m depth range) Richardson number.
An alternate derivation of $\sigma_{vi}$ is:
$ R_i = {N^2 \over { ( \frac {\partial u}{\partial z} )^2 + ( \frac {\partial v}{\partial z} )^2} }$
$ N^2 = {{-g \over \rho_0 }\frac {\partial \rho}{\partial z}} $
After hydrostatic and geostrophic approximations
$f \frac {\partial v}{\partial z} = {{-g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad f \frac {\partial u}{\partial z} = {{g \over \rho_0 }\frac {\partial \rho}{\partial y}} $
so
$\frac {\partial v}{\partial z} = {{-1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad \frac {\partial u}{\partial z} = {{1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial y}} $
$\therefore$
$ R_i = {f^2 N^2 \over \underbrace{{ {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial y})^2 + {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial z})^2} }_\text{m^4}}\quad\quad = \quad\quad {f^2N^2 \over m^4} $
$\sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}$
$RX_1 = RX_{east} = \Delta\rho_x = \rho_{i+1,j} - \rho_{i,j}$
$RY_1 = RY_{north} = \Delta\rho_y = \rho_{i,j+1} - \rho_{i,j}$
$RZ_1 = RZ_{k+1} = \Delta\rho_z = \rho_{k} - \rho_{k+1}$
$\displaystyle{1 \over L_{R_i}} \displaystyle\int_{2000m}^{100m} \left\lbrace { {-g\over\rho_0}{\frac {\partial \rho} {\partial z}} \over { {g^2 \over \rho_0^2 } \left[( \frac{\partial \rho}{\partial y})^2 + ( \frac{\partial \rho}{\partial z})^2\right] } } \right\rbrace dz$
Note: missing $f^2$ which will be cancelled when forming $\sigma_{vi}$
$\quad\quad$ so $\cdots$ this is not $R_i$
Implementation notes
For Level K
Numerator : Top $= -grav * RZ_{SAVE}(\cdots k+1) * dzwr(k)$
Denominator :
$ \begin{align} work1 = p25 & * ( RX(..,i_{east},k)^2 \\ & + RX(..,i_{west},k)^2 \\ & + RX(..,i_{east},k+1)^2 \\ & + RX(..,i_{west},k+1)^2 ) / DXT(i,j)^2 \\ \end{align} $
$ \begin{align} work2 = p25 & * ( RY(..,j_{north},k)^2 \\ & + RY(..,j_{south},k)^2 \\ & + RY(..,j_{north},k+1)^2 \\ & + RY(..,j_{south},k+1)^2 \\ \end{align} $
$work3 = {\left( TOP \over (grav^2*(work1+work2))\right)}*dzw(k)$
Notes:
1)Need to be careful at top and bottom of ocean
2)Accurate dzw(k) for each (i,j) to form $L_{R_i}$
3) When constructing $sigma$ itself, use $RZ_{SAVE}$ with a minimum N value
4) use eps2
# <!-- collapse=True -->
sigma_old = nc.variables['SIGMA_AVG'][0,:,:]
sigma_old_nearest = pyresample.kd_tree.resample_nearest(orig_def, sigma_old, \
targ_def, radius_of_influence=500000, fill_value=None)
sigma_avg1 = nc.variables['SIGMA_AVG1'][0,:,:]
sigma_avg1_nearest = pyresample.kd_tree.resample_nearest(orig_def, sigma_avg1, \
targ_def, radius_of_influence=500000, fill_value=None)
clevs=arange(.005,.05,.0025)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,sigma_avg1_nearest*86400,
levels1=clevs,
figsize=(14,5),
setunder='darkblue',
setover='darkred',
title="New Eady inverse time scale calc using alt sigma $(days^{-1}) \sim {{m^2} \over N}$")
yticks=[.01,.02,.03,.04,.05]
cbar.set_ticks(yticks)
tulloch_eady=Image(filename=basedir+'/steering/tulloch_eady.png')
display(tulloch_eady)
# <!-- collapse=True -->
################# PLOT u_rms with Leddy #######################
cm2m=.01
alpha=5.
u_rms_leddy = alpha*sigma_avg1_nearest*l_eddy_calc # units = 1*1/s*cm = cm/s
u_rms_leddy_mps=u_rms_leddy*cm2m
clevs=np.linspace(0,.4,40)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,u_rms_leddy_mps,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="$u_{rms} ({m \over s}) \quad alpha * \sigma * l_{eddy}$ (alpha=%i,max = %i/min = %i) "%(alpha,u_rms_leddy_mps.max(),u_rms_leddy_mps.min()))
#customize colorbar ticks
yticks=[0,.1,.2,.3,.4]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.2])
pylab.xticks([0,.1,.2])
# <!-- collapse=True -->
bates_urms=Image(filename=basedir+'/steering/bates2013-urms.png')
display(bates_urms)
# <!-- collapse=True -->
################# PLOT u_rms^2 convert m^2/s^2 #######################
u_rms_leddy_mps_sq=u_rms_leddy_mps*u_rms_leddy_mps
clevs=np.linspace(0,.06,40)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,u_rms_leddy_mps_sq,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="$u_{rms}^2 ({m^2 \over s^2})$ (alpha=%i,max = %g/min = %g) "%(alpha,u_rms_leddy_mps_sq.max(),u_rms_leddy_mps_sq.min()))
#customize colorbar ticks
yticks=[0,.02,.04,.06]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06])
# <!-- collapse=True -->
urms_bates=Image(filename=basedir+'/steering/Bates_urms_sq.jpg')
display(urms_bates)
# <!-- collapse=True -->
uminusc_sq_ovr_urms_sq = uminusc_mps_sq/u_rms_leddy_mps_sq #units cm/s / cm/s
uminusc_sq_ovr_urms_sq_masked=np.ma.masked_where(uminusc_sq_ovr_urms_sq >100., uminusc_sq_ovr_urms_sq, copy=True)
clevs=np.logspace(-3,4,21)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_sq_masked,
addzonal=True,
levels1=clevs,
norm=norm,
figsize=(14,5),
setover='darkred',
setunder='darkblue',
title="${(U-c)}^2 \over u_{rms}^2$ (masking values over 100)")
#customize colorbar ticks
yticks=np.logspace(-3,4,15)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,100])
pylab.xticks([0,50,100])
clevs=arange(0,5,.1)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_sq_masked,
addzonal=True,
levels1=clevs,
figsize=(14,4),
setover='darkred',
setunder='darkblue',
title="${(U-c)}^2 \over u_{rms}^2$ (masking values over 100)")
#customize colorbar ticks
yticks=[0,2,4]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,8])
pylab.xticks([0,4,8])
# <!-- collapse=True -->
bates_uminusc_sq_ovr_urms_sq =Image(filename=basedir+'/steering/Bates_suppress.jpg')
display(bates_uminusc_sq_ovr_urms_sq )
Parameterizing the Mixing term $\color{red}{L_{mix}}$
$\color{red}{L_{mix}} = \Gamma * L_{eddy} * Suppression$
$\Gamma = 0.35$
$b1 = 4$
$Suppression= {1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2 )}$
$\color{red}{L_{mix}} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)}$
# <!-- collapse=True -->
gamma=.35
b1=4.0
supp=1./(1+b1*uminusc_sq_ovr_urms_sq)
supp1b1=1./(1+1*uminusc_sq_ovr_urms_sq)
cm2m=.01
l_eddy_calc_m=l_eddy_calc*cm2m
lmix=gamma*l_eddy_calc_m*supp
lmix1b1=gamma*l_eddy_calc_m*supp1b1
clevs=np.logspace(0,5,40)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,lmix,
addzonal=True,
norm=norm,
levels1=clevs,
figsize=(14,5),
title="$L_{mix} = {\Gamma * L_{eddy} * Suppression} $")
pylab.xlim([0,10e3])
cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,3e4,1e5])
clevs=np.linspace(0,1,30)
norm=matplotlib.colors.Normalize()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,supp,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="Suppression factor ${1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2)}$")
#customize colorbar ticks
yticks=np.linspace(0, 1, 6)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.8])
pylab.xticks([0,.25,.5,.75])
# <!-- collapse=True -->
bates_suppression =Image(filename=basedir+'/steering/bates2013_suppression.png')
display(bates_suppression )
Parameterizing Eddy diffusivity $\color{red}{K}$
$\color{red}{K}=u_{rms}∗L_{mix}$
(b1, the suppression scaling factor, defaults to 4 but I also included a $\color{red}{K}$ plot with b1=1)
# <!-- collapse=True -->
kdiff=u_rms_leddy_mps*lmix
kdiff1b1=u_rms_leddy_mps*lmix1b1
clevs=np.logspace(0,4,40)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,kdiff,
addzonal=True,
norm=norm,
levels1=clevs,
figsize=(14,5),
title="Eddy diffusivity $K$ (log scale)")
#customize colorbar ticks
yticks=np.logspace(0, 4,5)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,11e3])
pylab.xticks([0,5e3,10e3])
clevs=np.linspace(0,1.e4,40)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,kdiff,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="Eddy diffusivity $K$ ")
#customize colorbar ticks
#cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,])
yticks=np.linspace(0, 10000,6)
print yticks
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,11e3])
pylab.xticks([0,5e3,10e3])
clevs=np.linspace(0,1.e4,40)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,kdiff1b1,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="Eddy diffusivity $K$ With a suppression scaling factor $b1$ of 1 instead of 4")
#customize colorbar ticks
#cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,])
yticks=np.linspace(0, 10000,6)
print yticks
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,11e3])
pylab.xticks([0,5e3,10e3])
# <!-- collapse=True -->
bates_k =Image(filename=basedir+'/steering/bates2013_K.png')
display(bates_k )
Steering Level status
1. L_eddy values ok
2. Rossby wave speed ok
3. Eady inverse timescale ok
a. Use alternate calculation otherwise has problems at equator and
b.Richardson averages are too high causing sigma to drop too low
4. Zonal Eddy Phase Speed (c)
a. Too high 10S:10N
b. Missing westward velocities in ACC which appear
in Hughes Data used in this parameterization
5. (U-c) Because of problem with c this term is also too high 10S:10N.
6. U_rms (alpha*sigma*l_eddy) OK
a. Using alpha of 5
7. Zonal mean velocity looks ok compared to ECCO annual average
8. Suppression OK
9. K is low by a factor of 5
Todo
1. Need to work on Zonal phase speed c
2. Plot vertical distribution of K
# <!-- collapse=True -->
urms_bates2=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-002.jpg')
urms_bates3=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-003.jpg')
urms_bates4=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-004.jpg')
urms_bates5=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-005.jpg')
urms_bates6=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-006.jpg')
urms_bates7=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-007.jpg')
urms_bates8=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-008.jpg')
urms_bates9=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-009.jpg')
urms_bates10=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-010.jpg')
urms_bates11=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-011.jpg')
display(urms_bates3,urms_bates4,urms_bates5,urms_bates6,urms_bates7,urms_bates8,urms_bates9,urms_bates10,urms_bates11)
CESM $U_{mean} = U_{resolved}$ integraged over {10,50,100,200,500} meters
# <!-- collapse=True -->
umean_cesm=Image(filename=basedir+'/steering/umean.jpg')
display(umean_cesm)
# <!-- collapse=True -->
bates_len_vel=Image(filename=basedir+'/steering/bates_len_vel.jpg')
display(bates_len_vel)
